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1.  Introduction 


Smart  projectiles  promise  to  provide  drastic  improvements  in  impact  point  perfonnance  of  artillery 
shells,  mortars,  grenades,  and  direct  fire  tank  rounds.  Although  the  potential  benefit  of  these  new 
weapons  is  quite  substantial,  the  technological  challenges  are  also  difficult.  In  contrast  to  un¬ 
guided  projectiles,  smart  projectiles  must  include  hardware  in  the  form  of  a  sensor  suite,  a 
microprocessor,  and  a  physical  control  mechanism  along  with  associated  software.  These  added 
components  must  be  small  to  fit  within  geometric  constraints  and  rugged  to  withstand  the  very 
harsh  launch  environment.  To  constrain  cost,  the  number  of  on-board  sensors  must  also  be 
minimized.  A  key  aspect  of  the  flight  control  system  software  is  the  logic  that  blends  available 
sensor  data  together  to  form  a  usable  feedback  signal  for  active  control.  This  is  the  classic  state 
estimator  or  observer  problem.  The  state  estimator  is  often  a  critical  element  within  a  smart 
weapon  flight  control  system  and  generally  affects  impact  performance  of  the  weapon  in  a  major 
way.  The  use  of  on-board  sensors  for  projectile  state  observation  is  plagued  by  the  noise  inherent 
in  physical  sensors,  while  the  use  of  on-board  state  propagators  is  plagued  by  inaccuracies  inherent 
in  simple  dynamic  models  (I). 

For  almost  50  years,  the  work  horse  of  state  estimation  has  been  the  celebrated  Kalman  filter  (2 
through  9),  which  allows  for  a  combination  of  imperfect  sensor  information  and  a  linear  model  of 
the  overall  system  to  be  blended  in  a  manner  that  provides  an  optimal  state  estimate  for  linear 
systems.  The  extended  Kalman  filter  (10,  11)  is  the  typical  way  that  estimators  are  designed  for 
nonlinear  physical  systems.  In  the  extended  Kalman  filter,  the  nonlinear  system  is  effectively 
approximated  by  a  linear  system,  and  the  standard  Kalman  filter  is  used  on  the  linear  system. 
Although  Kalman  filters  and  extended  Kalman  filters  are  used  extensively  in  practice,  they  are  not 
without  problems.  System  errors  and  measurement  errors  in  a  Kalman  filter  are  usually  assumed 
to  be  zero-mean  Gaussian  process  with  a  known  covariance.  Modeling  errors  seldom  assume  this 
form,  and  measurement  error  statistics  are  often  not  well  known,  leaving  users  to  select  error 
covariance  in  a  trial-and-error  fashion.  This  filter-tuning  process  can  be  time  consuming.  In 
some  applications,  it  can  be  difficult  to  get  the  estimates  to  converge. 

An  alternative  to  a  Kalman  filter  is  a  nonlinear  observer  that  works  directly  with  the  nonlinear 
dynamic  system  and  not  a  locally  linearized  approximation  of  the  system.  Because  of  a  large 
number  of  applications,  a  significant  literature  base  has  amassed  for  nonlinear  observers  (12,  IS, 
14).  Algebraic  methods  have  been  employed  to  transform  certain  classes  of  nonlinear  systems  into 
observer  normal  form,  where  observers  can  be  constructed  with  guaranteed  convergence 
properties.  Another  approach  to  construct  nonlinear  observers  lies  in  variable  structure  systems 
theory.  With  the  Lyapunov  theory,  extensions  of  standard  linear  observer  theory  can  be  made  to 
prove  asymptotic  stability  of  nonlinear  observers.  It  is  also  possible  to  employ  model  predictive 
control  techniques  to  form  nonlinear  model  predictive  observers.  The  focus  of  the  current  report  is 
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the  use  of  a  nonlinear  model  predictive  observer  to  estimate  the  state  of  a  smart  projectile.  The 
observer  is  based  on  a  two-step  process  consisting  of  an  initial  state  predictor  followed  by  a  state 
estimate  corrector.  We  generate  the  predictor  state  by  simply  integrating  the  projectile’s  fixed 
plane  equations  of  motion  forward  in  time  while  the  corrector  step  updates  the  state  estimate  by  the 
minimization  of  the  difference  between  measurements  and  model  predicted  measurements  over  a 
finite  duration  of  time  in  the  past.  Since  the  proposed  method  permits  a  general  nonlinear  repre¬ 
sentation  of  the  projectile,  the  corrector  step  requires  the  solution  of  a  nonlinear  minimization 
problem,  which  is  subsequently  solved  with  a  damped  Newton  procedure.  The  nonlinear  model 
predictive  observer  is  exercised  on  a  smart  projectile  that  employs  a  typical  array  of  sensors 
including  global  positioning  system  (GPS),  a  three-axis  gyroscope,  three-axis  accelerometers, 
and  a  three-axis  magnetometer.  Trade-off  studies  are  shown  that  vary  the  estimation  horizon. 


2.  Projectile  Flight  Dynamic  Model 


A  six-degree-of-freedom  (6-DOF)  rigid  projectile  model  has  been  proved  to  accurately  simulate 
the  dynamics  of  a  projectile  in  atmospheric  flight  (75).  The  6  DOF  are  comprised  of  the  three 
translational  components  describing  the  position  of  the  projectile’s  center  of  mass  and  the  three 
Euler  angles  describing  the  orientation  of  the  projectile  with  respect  to  a  fixed  inertial  axis. 
Figures  1  and  2  provide  a  schematic  of  the  degrees  of  freedom. 


Figure  1.  Position  definitions. 
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Figure  2.  Orientation  definitions. 


Although  commonly  expressed  and  computed  in  the  projectile  body  frame,  the  equations  of 
motion  for  the  6-DOF  model  are  expressed  in  the  no-roll  frame  to  reduce  computational 
requirements  within  the  observer. 
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The  force  acting  on  the  projectile  in  equation  3  is  comprised  of  the  weight  force  (W)  and  the 
aerodynamic  force.  The  aerodynamic  force  is  split  into  a  standard  (A)  and  Magnus  (M) 
aerodynamic  force.  The  combination  of  forces  is  expressed  in  equation  5. 
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Equation  6  provides  the  expression  for  the  weight  force  in  the  no-roll  coordinate  system. 
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Equation  7  provides  the  expression  for  the  aerodynamic  force  in  the  no-roll  coordinate  system. 
This  force  acts  upon  the  projectile  at  the  aerodynamic  center  of  pressure. 
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Equation  8  provides  the  expression  for  the  Magnus  force  in  the  no-roll  coordinate  system.  The 
Magnus  force  acts  upon  the  projectile  at  the  Magnus  force  center  of  pressure. 
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Equations  7  and  8  are  based  on  Mach  number-dependent  coefficients,  the  aerodynamic  angles  of 
attack,  and  the  total  aerodynamic  velocity. 


The  moment  acting  on  the  projectile  in  equation  4  is  comprised  of  the  moment  attributable  to  the 
standard  aerodynamic  force  (A),  the  moment  attributable  to  the  Magnus  aerodynamic  force  (M) 
and  the  unsteady  aerodynamic  moments  (ETA)  as  shown  in  equation  9. 
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The  moment  attributable  to  the  aerodynamic  force  is  expressed  in  equation  10. 
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The  moment  attributable  to  the  Magnus  force  is  expressed  in  equation  1 1 . 
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The  unsteady  aerodynamic  moments  acting  on  the  projectile  are  expressed  in  equation  12. 
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These  coefficients  are  projectile  specific  functions  of  the  Mach  number  of  the  projectile. 


3.  Sensor  Model 


Within  the  corrector  part  of  the  observer,  estimated  motion  of  the  projectile  is  contrasted  with  a 
given  sensor  suite  to  iteratively  determine  state  estimates.  Sensor  suites  on  a  smart  weapon 
typically  include  gyroscope,  accelerometer,  magnetometer,  and  GPS  measurements.  To  properly 
blend  these  data  together  in  the  observer,  an  analytical  sensor  model  is  required  for  each  sensor 
type,  which  relates  sensor  readings  to  projectile  motion  via  rigid  body  kinematics.  Major  error 
sources  such  as  noise  error,  bias  error,  scale  factor  error,  cross-axis  sensitivity/misalignment,  and 
sensor  misposition  can  significantly  affect  sensor  readings  and  must  therefore  be  included  in  the 
sensor  model.  The  acceleration  experienced  by  the  location  of  an  accelerometer  is 
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where 
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Also,  the  angular  velocity  of  the  gyroscopes  is 


SGx 

cy 

CG.v 

A 

P 

COy 

>  =  < 

>  +  < 

>  + 

cr 

Gy 

Sr 

Gy 

CGy 

[M 

q 

(14) 

°E . 

Pm. 

Pm. 

cx 

_CGz 

cy 

CGz 

SGz  J 

r 

Magnetometers  measure  the  dot  product  between  the  earth’s  magnetic  field  vector  and  the  sensi¬ 
tive  axis  of  the  magnetometer.  Although  they  do  not  directly  measure  any  rigid  projectile  model 
states,  three  magnetometers  can  be  arranged  orthogonally  so  that  they  project  the  earth’s  magnetic 
field  vector  onto  the  projectile  body  frame  and  aid  in  estimating  projectile  orientation. 
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GPS  measures  the  position  and  velocity  of  the  receiver  relative  to  the  earth,  which  is  given  by  the 
following: 
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4.  Nonlinear  Model  Predictive  Observer 


Consider  the  projectile  dynamic  system  just  described.  Using  the  state  vector  defined  by  x  =  {x,y, 
z(j),0,  i//,u,v,w,p,q,r  }T,  we  can  describe  the  dynamic  system  by  the  following: 

x  =  f(x)  (17) 

Assuming  that  several  quantities  are  measured  that  are  related  to  the  state  so  that  the  output  is 
described  by  the  following: 

y  =  h(x )  (18) 
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It  is  desired  to  estimate  the  state  x  at  any  time  with  the  use  of  the  math  model  and  measurements. 
The  state  estimate  will  be  called  x  and  the  measurement  reconstruction  uses  the  estimates 
v  =  h  (x) .  In  order  to  make  a  state  estimate,  a  quadratic  cost  function  is  minimized  over  a  discrete 

horizon  (see  figure  3). 
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where  W  is  positive  definite  and  k  is  the  selected  measurement  horizon. 

In  solving  the  minimization  problem,  it  is  desirable  to  express  the  cost  function  as  a  function  of 
the  state  estimate  at  j.  Equation  20  shows  the  cost  function  as  it  relates  to  the  current  state 
estimate. 
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g.  (xj  j  is  a  function  that  relates  the  state  estimate  at  the  current  time  to  past  sensor  measurements. 
Its  structure  is  discussed  later  in  this  report. 


For  the  state  estimate  to  be  an  optimal  solution,  the  cost  function  must  satisfy  the  gradient  condition. 
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In  order  to  satisfy  equation  21,  a  numerical  technique  is  used  because  J  is  highly  nonlinear.  A 
Newton  type  method  is  used  to  solve  for  x  that  will  drive  dJ/dx  =  0  .  To  achieve  this  condition 
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numerically,  J  is  approximated  as  a  local  quadratic  function  (equation  22),  and  a  perturbation 
on  the  state  is  solved  for  so  that  the  gradient  of  the  local  model  is  zero  (equation  23). 
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The  line  search  parameter  a  is  implemented  in  order  to  accelerate  convergence  (equation  24): 
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When  a  set  of  nonlinear  equations  is  solved  iteratively,  a  starting  guess  is  required.  This  initial 
guess  is  generated  by  integrating  the  equations  of  motion  (1  through  4)  forward  one  time  step. 
This  can  be  viewed  as  the  predictor  step,  and  solving  the  nonlinear  equations  can  be  viewed  as 
the  corrector  step. 

A  key  aspect  of  this  technique  is  linking  the  output  at  one  point  in  time  to  the  state  at  the  time  of 
interest.  To  create  this  link,  the  dynamic  equations  of  motion  can  be  used  and  applied  with  a 
suitable  discretization.  For  this  problem,  Euler’s  formula  was  used  for  its  simplicity  and  ease  of 
computation. 

*,-i  =Xj  -Atf(xj)  (25) 

This  link  can  be  used  recursively  to  link  the  state  at  any  two  times. 

Xj-2  =  Xj-1  -  Atf  (xj- 1  )  =  Xj-  Atf  (xj )  -  At/  (xj  -  A tf  (xj )) 

Xj.  3  =  Xj  -  A  tf  (5tj )  -  A  tf  (xj  -  A  tf  (x. ))  -  A  tf  (x.  -  A  tf  (xy )  -  A  tf  (x;  -  A  tf  (/.))) 


Equation  24  requires  the  partial  derivatives  of  the  measurements.  These  derivatives  can  be  com¬ 
puted  numerically,  partially  numerically,  and  partially  analytically,  or  completely  analytically. 
Using  the  analytical  approach,  we  get  the  following  equations: 


dJ  k+l 
fix.. 


fiv 

fix; 


z 

(=1 

k+l 

=z 

1=1 


dy(xt)  T 


fix,. 


V'(y(xl)-y(x,))  +  (y(x,)-y(x,)jr  W 


) 


fiX; 


d2y(Xi) 


dx2 


w{y{xi)-y{xi))  +  2 


fiy(T;.) 


fix,. 


(27) 


dxj 


fix ; 


The  second  partial  derivatives  used  in  equation  27  can  typically  be  neglected  because  they  are  very 
small  when  compared  with  the  first  partial  derivative  term  (16). 
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Equation  27  requires  more  derivatives  to  be  determined. 


jfei) 

dxj 


dy(*J- 1) 

fix,. 


dh 

fix 

dh 

dx 


X  <9x, 
xj- 1  J 


(xj-Atf^j)) 


f  ^ 

I-At^ 


fix 


v 


These  derivatives  can  also  be  recursive. 

dv(*;-2)  _  ^(^-2)  dXj_ 2  fix,.,, 

fix,.  fix, ,  fix,,  fix, 

J  J  A  J  1  J 


(*/-*) 


fix. 


dh 

fix 


r 


n 

ij-N  P= 1 


/-At 


fix 


lj-N+p  J 


(28) 


(29) 

(30) 


5.  Sample  Results 


In  order  to  exercise  the  methodology  developed  here,  simulation  results  are  generated  for  a  sample 
system.  The  projectile  considered  is  a  typical  155-mrn  spin-stabilized  projectile.  The  projectile  is 
assumed  to  have  a  sensor  suite  on  board  which  includes  GPS,  a  three-axis  accelerometer,  a  three- 
axis  magnetometer,  and  a  three-axis  gyroscope. 

Convergence  of  the  corrector  portion  of  the  observer  is  shown  in  figure  4  for  a  single  instant  in 
time  for  a  single  magnetometer  axis.  In  figure  4,  a  total  of  15  points  is  used  in  the  estimation 
horizon,  and  the  numerical  solution  exhibits  rapid  convergence.  The  cost  J  with  a  guess  at  initial 
conditions  was  7.1735e+004  and  was  reduced  to  1.9027977e-001  and  1.9027976e-001  after 
estimates  4  and  5,  respectively.  The  cost  had  not  only  converged,  but  the  estimated  output  was 
very  acceptable  for  the  entire  horizon.  These  results  are  typical  for  all  measurements. 

A  trade  study  was  conducted  to  determine  an  appropriate  number  of  points  in  the  measurement 
horizon  k.  A  number  of  various  horizons  were  selected,  ranging  from  three  data  points  to  42  data 
points.  The  measurement  data  were  populated  with  noise  with  standard  deviations  shown  in  table  1 
and  remained  the  same  for  each  case  in  the  study.  Also,  the  initial  conditions  at  j  =  1  were  the  same 
for  each  case  of  the  study.  The  results  of  this  trade  study  are  shown  in  the  figure  5.  As  the  number 
of  points  in  the  prediction  increases,  the  sensor  noise  becomes  increasingly  filtered.  From  this  trade 
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study,  it  was  determined  that  a  horizon  of  2 1  data  points  gave  satisfactory  estimates  without 
requiring  intense  computations. 


Figure  4.  Single  time  step  iterations. 

Table  1.  Sensor  noise  standard  deviation. 


Sensor 

Noise  Standard 
Deviation 

Sensor 

Noise  Standard 
Deviation 

GPS 

Accelerometer 

X 

4.8  ft 

X  -axis 

0.03  ft/sec2 

y 

4.8  ft 

y  -axis 

0.03  ft/sec2 

z 

10.8  ft 

z  -axis 

0.03  ft/sec2 

X 

0.6  ft/sec 

Magnetometer 

y 

0.6  ft/sec 

X  -axis 

0.001 

z 

0.6  ft/sec 

y  -axis 

0.001 

Gyroscope 

z  -axis 

0.001 

P 

0.1  rad/sec 

q 

0.1  rad/sec 

r 

0.1  rad/sec 
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Figure  5.  Trade  study  results. 

To  further  illustrate  how  greatly  this  estimator  can  reduce  to  amplitude  and  frequency  of  noise, 
figure  6  compares  the  sensor  output,  estimated  state,  and  actual  state  for  a  single  GPS  signal. 

With  a  measurement  horizon  of  21  points,  state  estimates  were  found  over  a  full  projectile 
trajectory.  The  comparisons  of  each  of  the  states  with  the  state  estimates  are  plotted  in  figures  7 
through  22.  Also  plotted  are  the  some  of  the  errors  between  the  actual  state  and  the  state  estimate 
along  the  trajectory. 


6.  Conclusions 


A  nonlinear  model  predictive  observer  specialized  to  smart  projectiles  reported  here  can  success¬ 
fully  predict  the  full  state  vector  of  a  projectile  trajectory.  The  observer  is  based  on  a  two-step 
process  consisting  of  an  initial  state  predictor  followed  by  a  state  estimate  corrector.  We  generate 
the  predictor  state  by  simply  integrating  the  projectile’s  fixed  plane  equations  of  motion  forward  in 
time  while  the  corrector  step  updates  the  state  estimate  by  the  minimization  of  the  difference 
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between  measurements  and  model  predicted  measurements  over  a  finite  duration  of  time  in  the  past. 
Since  the  method  permits  a  general  nonlinear  representation  of  the  projectile,  the  corrector  step 
requires  the  solution  of  a  nonlinear  minimization  problem  which  is  subsequently  solved  with  a 
damped  Newton  procedure.  The  nonlinear  model  predictive  observer  is  exercised  on  a  smart 
projectile  that  employs  a  typical  array  of  sensors  including  GPS,  a  three-axis  gyroscope,  three-axis 
accelerometers,  and  a  three-axis  magnetometer.  Results  with  this  smart  weapon  observer  are 
promising  since  the  observer  is  capable  of  estimating  the  full  state  of  the  projectile,  including 
orientation  and  translational  velocity,  when  realistic  noise  and  bias  errors  are  included  in  the  sensor 
measurements.  Relative  to  other  existing  observers,  however,  the  proposed  observer  is  computa¬ 
tionally  intensive  because  of  the  need  to  solve  a  nonlinear  minimization  problem  at  each  compu¬ 
tation  cycle  of  the  observer. 


Figure  6.  Noise  filtering. 
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Figure  7.  Projectile  range. 
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Figure  8.  Projectile  cross  range. 
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Figure  9.  Projectile  altitude. 
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Figure  10.  Projectile  roll  angle. 
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Figure  11.  Projectile  Euler  angle  9. 


Figure  12.  Projectile  Euler  angle  \|/. 
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Figure  13.  Projectile  velocity  u. 
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Figure  14.  Projectile  velocity  v. 
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Figure  15.  Projectile  velocity  w. 


Figure  16.  Projectile  rate  p. 
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Figure  17.  Projectile  rate  q. 
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Figure  18.  Projectile  rate  r. 
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Figure  19.  Estimate  error  in  range. 
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Figure  20.  Estimate  error  in  <j>. 
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Figure  21.  Estimate  error  in  u. 


Figure  22.  Estimate  error  in  p. 
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Nomenclature 


x,y,z 

<t>,o,y/ 

ii,  v,  w 
P,q,r 
I 


aNi 


" Gi 


Ni 


G, 

M 

M 

R. 


Ni 


'  ©— 
SAi 


’Gi 


a 


B/I 


5JSi 


co 


B!I 


coB 


CO 


Ni 


(Oi 


Inertial  positions  of  the  system  mass  center. 

Euler  roll,  pitch  and  yaw  angles. 

Velocity  components  of  mass  center  in  a  no-roll  reference  frame. 
Angular  velocity  components  in  a  no-roll  reference  frame. 

Inertia  matrix  of  system. 

Bias  of  ith  accelerometer. 

Noise  of  ith  accelerometer. 

Acceleration  reading  of  ith  accelerometer. 

Cross  axis  sensitivity  of  ith  accelerometer  to  the  jth  direction. 
Cross  axis  sensitivity  of  ith  gyroscope  to  the  jth  direction, 
ith  component  of  the  Earth’s  magnetic  field. 

Bias  of  ith  GPS. 

Noise  of  ith  GPS. 

Reading  of  ith  GPS. 

Noise  of  ith  magnetometer. 

Reading  of  ith  magnetometer. 

Rotation  matrix  from  body  frame  to  sensor  frame. 

Distance  vector  from  projectile  mass  center  to  ith  accelerometer. 
Scale  factor  of  ith  accelerometer. 

Scale  factor  of  ith  gyroscope. 

Angular  acceleration  of  projectile  with  respect  to  the  ground. 
Misposition  of  ith  accelerometer  in  the  jth  direction. 

Angular  velocity  of  projectile  with  respect  to  the  ground. 

Bias  of  ith  gyroscope. 

Noise  of  ith  gyroscope. 

Angular  velocity  reading  of  ith  gyroscope. 

Vector  component  extraction  operator  for  frame  B  . 

Skew  symmetric  cross  product  operator  for  frame  B  . 
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